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' Abstract 

I A comparative study of simulated air shower longitudinal profiles is presented. An 

J> ■ appropriate thinning level for the calculations is first determined empirically. High 

statistics results are then provided, over a wide energy range, (lO 14,0 to 10 20 5 eV), 
for proton & iron primaries, using four combinations of the MOCCA & CORSIKA 
program frameworks, and the SIBYLL & QGSjet high energy hadronic interaction 
models. These results are compared to existing experimental data. The way in which 
the first interaction controls X max is investigated, as is the distribution of X max . 
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1 Introduction 



The energy spectrum of cosmic rays is a power law with the flux falling by 
three orders of magnitude for each decade increase in energy. At ~ 10 14 eV the 
flux becomes so low that current balloon and satellite experiments lack the 
exposure required to detect a significant sample of events. This is unfortunate 
as the nature of the primaries remains of great astrophysical interest. Where 
direct measurements are possible the cosmic rays are known to be mostly 
protons and atomic nuclei. The most plausible acceleration site is at the shock 
fronts produced by supernova explosions. However, theoretical considerations 
predict a maximum energy from this process of ~ 10 15 eV, whereas the energy 
spectrum is observed to continue with only small deviations up to > 10 19 eV. 
The origin of the particles at > 10 15 eV is somewhat mysterious. 

It has long been supposed that insight would result if the composition of the 
primaries could be measured. Due to the extremely low flux the only way to get 
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information on these particles is to study the extensive atmospheric cascades 
which they initiate. When a cosmic ray enters the atmosphere it collides with 
the nucleus of an air atom, producing a number of secondaries. These go on 
to make further collisions, and the number of particles grows. Eventually the 
energy of the shower particles is degraded to the point where ionization losses 
dominate, and their number starts to decline. 

It is a coincidence that at the energy where direct detection of the cosmic rays 
becomes impractical, the resulting air showers become big enough to be easily 
detectable at ground level. The number of particles in the cascade also be- 
comes large enough that the longitudinal profile, or number of particles versus 
atmospheric depth, becomes a smooth curve, with a well defined maximum. 
This maximum depth, referred to as X max , is often regarded as the most basic 
parameter of an air shower, and much effort has been expended to measure 
and interpret it. 

The depth of maximum increases with primary energy as more cascade gen- 
erations are required to degrade the secondary particle energies. For given 
total energy X max is related to the energy per nucleon of the primary. To first 
order the interactions occur between individual nucleons of the primary, and 
the target air nuclei. Therefore a shower initiated by a compound nucleus can 
be thought of as the superposition of many proton initiated showers, with 
correspondingly lower energy. 

Unfortunately, of course, the detail is not so simple. For a number of reasons 
extracting information on the nature of the cosmic ray primaries from the air 
showers they produce has proved to be exceedingly difficult. The most fun- 
damental problem is that the initial interactions are subject to large inherent 
fluctuations. This limits the event-by-event mass resolution of even an ideal 
detector. However, progress can still be made by looking at mean parameter 
values, or better, their distributions. 

The second major problem is that sophisticated modeling is required to predict 
the absolute value of an observable parameter which is expected for a primary 
of given type and energy. Nucleus-nucleus interactions at the energies of the 
first few cascade steps are well beyond the reach of accelerator experiments. 
Therefore it is necessary to rely on hadronic interaction models which attempt 
to extrapolate from the available data using different mixtures of theory and 
phenomenology. The lower energy part of the cascade can be modeled using 
well known physics, although the programs are complex with corresponding 
scope for errors. 

The depth of shower maximum has been determined by a number of experi- 
ments. In the energy range 10 14 to 10 16 eV it has been measured with varying 
degrees of directness using Cerenkov light [1-3]. The range 10 17 to 10 19 eV 
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Program frame 


High E Hadronic 


Low E Hadronic 


Electromagnetic 


IVKJUL/A 


Internal 
SIBYLL 


Internal 


Internal 


CORSIKA 


SIBYLL 
QGSjet 


GEISHA 


EGS4 



Table 1 



Summary of the four program-frame / interaction-model combinations. 

has been observed rather directly by the Fly's Eye detector through fluores- 
cence light [4]. Finally the region above 10 19 eV is the focus of the HiRes Fly's 
Eye [5], Auger Project [6], and others. In the past the experimental resolution 
and statistics have often been so poor that the mean value of X max has been 
discussed rather than its distribution — this is changing. 

The simulations required to interpret the data from any given experiment have 
usually been performed only for the energy range accessible to it. This is unfor- 
tunate since checking the consistency of experiments in adjacent energy ranges 
is critically important, given the uncertainty of the high energy hadronic in- 
teraction models. Additionally the exact value of X max for a given model can 
depend on the way in which the longitudinal profiles are recorded and fit. The 
purpose of this paper is to provide X max values with good statistical preci- 
sion, over a wide energy range, and computed in a consistent way using several 
hadronic models and two different cascade "framework" programs; for a more 
detailed discussion see [7]. 

The process of air shower simulation can be broken up into several parts. A 
framework program is required which handles the mechanics of the process 
and calls appropriate subroutines to model the interaction and propagation 
of the particles. Some fraction of the required transport and interaction mod- 
eling may be provided using third party code. In this paper two air shower 
simulation packages which have been heavily used in the literature are consid- 
ered. The first is MOCCA written by Hillas [8]. This originally used a simple, 
built-in hadronic interaction model, but has also been linked to SIBYLL [9]; 
all other modeling is handled internally The second program is CORSIKA, a 
well documented and thorough program prepared originally for the Kascade 
experiment [10]. It is linked to a number of high energy hadronic models, two 
of which are suitable for use over the very wide energy range of this study; 
SIBYLL and QGSjet [11]. An attractive feature of this program is the use of 
the well established High Energy Physics codes EGS4 [12] and GEISHA [13], 
for the electromagnetic, and lower energy hadronic modeling respectively. See 
Table 1 for a summary. 

Due to the inherent limitations of air shower fluctuations, and also because of 
poor experimental resolution and statistics, X max data is often compared only 
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to simulated values for proton and iron nuclei primaries. These are generally 
regarded as the extreme ends of the possible range. At lower energies the 
composition of cosmic rays tracks the general abundances of solar system 
material, with some modifications due to propagation spallation effects. Iron 
is the heaviest significantly abundant element. 



2 Technical Details 

At the highest cosmic ray energies it is absolutely necessary to use techniques 
which accelerate the simulation process. A popular approach is called thinning: 
below a threshold energy only a sub-set of the particles are tracked, with 
weightings to compensate for those discarded [8]. The threshold is usually 
specified as a fraction of the primary energy, and referred to as the "thinning 
level". 

For this study MOCCA92 and CORSIKA 5.62 were used. This version of 
CORSIKA includes a very similar thinning algorithm to MOCCA. In all cases 
the high energy hadronic interaction model was used with the set of inelastic 
cross sections provided by its authors. Electromagnetic particle energy cutoff 
was a uniform 0.2 MeV. 

When considering gross quantities, like the depth of shower maximum, it is 
possible to run the simulation codes with very severe thinning, and still obtain 
results of adequate quality. This means that many showers can be generated, 
over a multidimensional grid of primary parameters and shower models, within 
an acceptable computing time. The thinning process leads to longitudinal 
profiles which have large non-statistical fluctuations. The magnitude of these 
fluctuations increases with the severity of the thinning; see Figure 1. 

To recover the depth of shower maximum from "noisy" thinned profiles it is 
customary to fit them to an empirical cascade shape function. This is also 
necessary when analyzing experimental data as the quality is often poor. The 
following function was introduced by Gaisser and Hillas [14] as a "simple 
analytic parameterization" of the longitudinal profile of air showers: 

(V \ Xmax I A 
ir -\ exp(X max -X)/X, (1) 

-A- max ' 

where X is the atmospheric depth (in g cm -2 ), N max is the number of cascade 
particles at shower maximum, X max is the depth of maximum, and A is a 
characteristic length parameter (in the above reference fixed at 70 g cm -2 ). 
This is a Gamma Function, a form which naturally arises in cascade theory, 
and assumes that the first interaction is at X = 0. 
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Fig. 1. The longitudinal profiles of two simulated showers. At left the thinning level 
is 1CP 5 ' 5 of primary energy, and at right 1CP 3 ' 5 . 



A fourth parameter X is often introduced into Equation 1, ostensibly to allow 
for a variable first interaction point, 

N(X) = N max °— exp (X max - X)/\. (2) 



This seems somewhat inelegant as varying X does not correspond to a transla- 
tion along the X axis, unless X max is also changed. The following (equivalent) 
form is physically clearer, where X rise is the distance from first interaction to 
shower maximum, 

N{X)=N max ( X ~ X ° ) exp(X„ se + X -X)/A. (3) 



In practice, when fitting simulated hadronic cascade profiles to either Equa- 
tion 2 or 3, it turns out that X correlates poorly with the actual depth of 
first interaction, and frequently takes unphysical negative values. Experiments 
were made performing the fit with Xq fixed at the actual depth of first interac- 
tion X\. This produces a significantly poorer goodness of fit, and reduces the 
X max results by < 10 g cm -2 . The choice of Equation 2 or 3 also influences the 
X max results. For the remainder of this paper, in the interests of compatibility 
with other published results, Equation 2 was used with all four parameters 
free. Note that X is best regarded as simply an additional arbitrary shape 
parameter. 

To determine an appropriate thinning level for this study, sets of 500 proton 
showers were generated and fit, at each of 5 thinning levels; 10~ 3,5 , 10 -4,0 , 
10~ 4 ' 5 , 10~ 5 and 10~ 5,5 . Taking the 10~ 5,5 thinned distributions as reference, 



a Kolmogorov testQ was performed on each of the more heavily thinned sets, 
and for each of the fit parameters. This is a statistical test of the compati- 
bility in shape between two histograms — it yields the probability that they 
are drawn from the same parent distribution. All the fit parameters returned 
a high probability at thinning levels of 10 -5 and 10~ 4 ' 5 , i.e. the results are 
indistinguishable within the statistics of a 500 event set. X max itself is extraor- 
dinarily robust, remaining unbiased even at 10~ 3 ' 5 thinning. To be conservative 
a value of 10~ 4,5 was selected for the main study. 

The compatibility of the parameter distributions was also checked when vary- 
ing the zenith angle of the primary from to 60 deg. X max was unaffected, 
but interestingly N max showed a small systematic increase with angle [7]. 



3 X max Results 

3.1 Mean value of X max 

For the main study sets of 500 events were run at 14 half decade energy steps 
between 10 14 and 10 20 ' 5 eV, with the 4 combinations of framework program 
and high energy hadronic interaction model given in Table 1. Sets were gener- 
ated with both proton and iron nucleus primaries. This gives 500 x 14 x 4 x 2 = 
56, 000 showers. The showers were run at a thinning level of 10~ 4 5 of primary 
energy, and a zenith angle of 45 deg. The resulting profiles were fit to Equa- 
tion 2. Figure 2 shows the mean value of X max plotted against energy over the 
complete range; numerical values are given in Table 3. MOCCA-SIBYLL and 
CORSIKA-SIBYLL proton results are in good agreement, and the iron results 
are also close. This is very encouraging — the framework programs are com- 
plex and entirely independent — nevertheless they produce the same result. At 
higher energies the older MOCCA-Internal model diverges strongly to deeper 
Xmax- CORSIKA-QGSjet produces much shallower results than SIBYLL at 
all energies; so much so that at 10 20 5 eV MOCCA-Internal iron is equal to 
CORSIKA-QGSjet proton. 

Figure 3 shows a comparison between published mean X max data from two ex- 
periments and the CORSIKA calculations. The data for E < 10 17 eV are from 
the BLANCA experiment [3], and for E > 10 17 eV from the Fly's Eye [4]Q. 
The Fly's Eye data contains a small un-corrected experimental bias, the re- 
moval of which would shift the lower energy points higher in the atmosphere 



1 The test as implemented in the HBOOK function HDIFF [15] was used. 

2 These points are used rather than the newer ones in [16] since they have a much 
wider energy range, and their quoted errors are comparable, or better. 
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Fig. 2. The mean depth of shower maximum versus energy for p^roion^and iron 
primaries, and four program-frame / interaction-model combinations. 

by around 20 g cm -2 [17]. Both SIBYLL and QGSjet are consistent with the 
data, in that the value remains within the proton-iron bounds. Also the gen- 
eral trends in composition versus energy are the same under the two models, 
although the absolute value and size of the changes differ. There is some evi- 
dence at ~ 10 16 eV that QGSjet is a more realistic model than SIBYLL [3,18]; 
the extrapolation to the highest energies in both models must be regarded as 
tentative at best. 



3.2 Influence of the First Interaction on X„ 



Why do the results shown in Figure 2 differ between the models? The proton- 
air cross sections used by SIBYLL and QGSjet are sufficiently similar that the 
mean free paths differ by < 5 g cm" 2 over the complete energy range. This is 
to be compared to the 30-50 g cm -2 difference in mean X max . 
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Fig. 3. A Comparison of experimental mean X„ 
high energy hadronic interaction models. 



data with simulations using two 



When investigating the way in which the first interaction controls X max it is 
natural to subtract out the position of the first interaction; X r i se = X max — X\. 
Elasticity is defined as the energy fraction of the most energetic secondary. 
Normally a large fraction of the primary energy continues in the form of a 
"leading nucleon" and the remainder is split between many secondary pions 
and nucleons. Figure 4 shows X T i se versus first interaction elasticity at 10 19 eV. 
For events where the first interaction is catastrophic (small elasticity), the re- 
sulting shower takes fewer generations to reach maximum, and the correlation 
is strong. As elasticity becomes larger, the first interaction is no longer the 
controlling factor, and the relationship weakens. Interestingly, both models 
exhibit approximately the same correlation between elasticity and X r i se . 



Figure 5 shows the elasticity distributions versus energy. The reason the mod- 
els produce different mean X max values is primarily because of their different 
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Fig. 4. The correlation of X r i se with first interaction elasticity for primary protons 
at 10 19 eV. 
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Fig. 5. Elasticity of proton-air interactions versus energy for two models. The side 
length of the boxes is proportional to the bin contents. 



elasticity distributions; QGSjet produces many more "hard" events, which lead 
to less deeply penetrating showers. SIBYLL has rather constant behaviour ver- 
sus energy, while QGSjet is a more radical model, showing a stronger change; 
this is why the corresponding mean X max results diverge with increasing en- 
ergy. 

Hadronic interactions models are complex and esoteric, with many parameters 
which can potentially be compared. The significance of Figure 4 is the realiza- 
tion that, for calculations of X max at least, the most important characteristic 
is also a very simple one: how much of the primary energy is expended in the 
first interaction? 
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Fig. 6. The distribution of X max versus energy shown as a band containing the 
central 68% of the distribution. 

3.3 Distribution of X max 

In Figure 2 it can be seen that QGSjet predicts that the difference between 
proton and iron mean X max decreases significantly with energy. However, the 
fluctuations do not decrease correspondingly, and the proton and iron dis- 
tributions overlap to an increasing extent. If this model is correct, greater 
experimental statistics would be required to determine the mean composition 
with given accuracy. The situation is illustrated in Figure 6 where the bands 
contain 68% of the data (i.e. spanning the 16% and 84% points of the integral 
distribution). CORSIKA-SIBYLL shows stronger separation improving with 
increasing energy, while CORSIKA-QGSjet has weaker separation degrading 
with energy. 

Proton primaries are deflected less by magnetic fields than more highly charged 
particles of the same total energy. It has been suggested that attempts to lo- 
cate the origin of the highest energy cosmic rays, by studying their arrival 
directions, could be enhanced by making cuts on composition sensitive pa- 
rameters, to increase the fraction of protons in the data sample. This would 
clearly work much less well if QGSjet is a more realistic model than SIBYLL. 

For proton primaries the distribution of X max is strongly asymmetric, with a 
tail to deep X max . This is presumably the result of fluctuations in the first 
interaction point, and is therefore connected to the proton- air cross section, 
which is a quantity of fundamental interest. Earlier simulations have suggested 
that, 

A c Ap_ a j r , (4) 
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Fig. 7. An example X max distribution with exponential trailing edge fit. 



Model 


Ratio c 


X 2 /ndf 


MOCCA-Internal 


1.32 ±0.03 


0.77 


MOCCA-SIBYLL 


1.16 ±0.03 


1.24 


CORSIKA-SIBYLL 


1.15 ±0.03 


1.38 


CORSIKA-QGSjet 


1.30 ±0.04 


0.43 



Table 2 



Decrement / mean-free-path data. 



where A is the exponential slope, or "decrement" , of the trailing edge of the 
X max distribution, X p - a i r is the proton-air mean free path, and c is a constant 
of proportionality with a value between 1.2 and 1.6 dependent on hadronic 
interaction model [19]. 

The X max distributions were fit to an exponential starting at 100 g cm -2 be- 
yond the peak. An example distribution, with the fit, is shown in Figure 7. To 
avoid biasing the results it is essential to use a maximum likelihood algorithm 
since the bins on the far tail necessarily contain few events. Figure 8 shows the 
value of the ratio c, plotted versus energy, for one of the models. Testing each 
set of results against the hypothesis of energy independence yields the values 
given in Table 3.3. With the available statistics, the reduced y 2 numbers show 
little evidence for energy dependence. The SIBYLL based models give values 
close to 1.15, while the other two are around 1.30; this difference appears to 
have significance. 



11 



o 



:ORSIKA-SIBYLL 



14 15 16 17 



19 20 

Log 1D (E (eV)) 



Fig. 8. Ratio c plotted versus energy. The most probable value is shown as a hori- 
zontal line. 

4 Conclusions 



When running shower simulations to study X max it is better to generate heav- 
ily thinned showers, with explicit low energy hadronic and electromagnetic 
cascades, than to rely on analytic approximations for these parts of the cal- 
culation. The latter has frequently been done in the past, leading to concerns 
that the results are biased to an unknown extent. When working with an 
explicit, but thinned, cascade simulation it is possible to determine an appro- 
priate thinning level empirically, by comparing against more lightly thinned 
results. 

Carefully calculated X max results have been presented, over a wide energy 
range, for proton & iron primaries, using four combinations of framework 
program and high energy hadronic interaction model. It is hoped that these 
will be of use for future comparisons with experimental data, and with other 
simulation results. 

The way in which the first interaction controls X max has been investigated. The 
influence is strong — if one were to use model A for the first few interactions, 
and model B thereafter, the mean X max results would be close to using model 
A throughout. QGSjet predicts that the separation between proton and iron 
X m ax declines at the highest energies. If this is true it is unfortunate from an 
experimental perspective. 

It would be very useful if a common reference set of showers were made avail- 
able by the authors of new, or modified, hadronic interaction models. For the 
purposes of longitudinal profile comparison the set used here seems adequate; 
the raw and processed output is available online [20]. 
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Table 3 ~ 

Mean and standard deviation values of X max for proton and iron primaries, and 
four program-frame / interaction-model combinations. Each pair of numbers comes 
from a 500 shower set. \a 



